packages <- c("CIMseq", "CIMseq.testing", "tidyverse", "ggthemes")
purrr::walk(packages, library, character.only = TRUE)
rm(packages)
##DATA
load('../data/CIMseqData.rda')
load('../data/sObj.rda')
if(!dir.exists('../figures')) dir.create('../figures')
p <- plotUnsupervisedClass(cObjSng, cObjMul, palette('total'))
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.classes.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)
p <- plotUnsupervisedMarkers(
cObjSng, cObjMul,
c("Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1"),
pal = RColorBrewer::brewer.pal(8, "Set1")
)
ggsave(
plot = p,
filename = '../figures/MGA.enge.markers.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)
p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Mki67")
ggsave(
plot = p,
filename = '../figures/MGA.enge.Mki67.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)
p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Hoxb13")
ggsave(
plot = p,
filename = '../figures/MGA.enge.Hoxb13.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)
classHeatmap(
data = data.frame(
gene = markers$gene,
class = markers$cluster
),
counts.log = getData(cObjSng, "counts.log"),
classes = tibble(
sample = colnames(getData(cObjSng, "counts")),
class = getData(cObjSng, "classification")
)
)
Cell type markers
markers <- c(
"Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1", "Mki67",
"Hoxb13", "Plet1", "Osr2", "Reg4", "Sval1"
)
gene.order <- c(
"Hoxb13", "Osr2", "Atoh1", "Reg4", "Sval1", "Plet1", "Mki67", "Lgr5", "Alpi",
"Slc26a3", "Lyz1", "Dclk1", "Chga", "Ptprc"
)
d <- getData(cObjSng, "counts.cpm")[markers, ] %>%
t() %>%
as.data.frame() %>%
rownames_to_column("Sample") %>%
as_tibble() %>%
gather(gene, cpm, -Sample) %>%
inner_join(tibble(
Sample = colnames(getData(cObjSng, "counts")),
Classification = getData(cObjSng, "classification")
)) %>%
mutate(Classification = parse_factor(Classification, levels = classOrder.MGA('total'))) %>%
mutate(gene = parse_factor(gene, levels = gene.order)) %>%
arrange(Classification, gene) %>%
mutate(Sample = parse_factor(Sample, levels = unique(Sample))) %>%
mutate(type = case_when(
str_detect(Classification, "^S") ~ "Small intestine",
str_detect(Classification, "^C") ~ "Colon",
TRUE ~ "Other"
)) %>%
mutate(type = parse_factor(
type,
levels = c("Small intestine", "Colon", "Other")
))
## Joining, by = "Sample"
under <- d %>%
mutate(idx = 1:nrow(.)) %>%
group_by(Classification) %>%
summarize(
min = (min(idx)), max = (max(idx) - 1), median = median(idx)
) %>%
ggplot() +
geom_segment(
aes(x = min, xend = max, y = 0, yend = 0, colour = Classification)
) +
geom_text(
aes(median, y = 0, label = Classification),
angle = 90, nudge_y = -0.01, hjust = 1, colour = "grey10"
) +
scale_colour_manual(values = palette('total')[classOrder.MGA('total')]) +
ylim(c(-1, 0)) +
theme_few() +
scale_x_continuous(expand = c(0, 0)) +
theme(
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
panel.border = element_blank()
) +
guides(colour = FALSE)
over <- d %>%
ggplot() +
geom_bar(aes(Sample, cpm, fill = Classification), stat = "identity") +
facet_grid(gene ~ type, scales = "free", space = "free_x") +
scale_fill_manual(values = palette('total')[classOrder.MGA('total')]) +
theme_few() +
theme(
axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
strip.background.x = element_rect(fill = "white")
) +
labs(y = "CPM") +
guides(fill = FALSE)
library(patchwork)
p <- over + under + plot_layout(ncol = 1, heights = c(2, 1))
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.geneBar.pdf',
device = cairo_pdf,
height = 300,
width = 300,
units = "mm"
)
Stemness dotplot
load('../data/seurat.rda')
markers <- c("Wnt3", "Wnt2b", "Wnt4", "Wnt5a", "Lgr5", "Axin2", "Ascl2")
gene.order <- markers
#scale.func <- scale_size
scale.func <- scale_radius
scale.min <- NA
scale.max <- NA
p <- getData(cObjSng, "counts.cpm")[markers, ] %>%
t() %>%
as.data.frame() %>%
rownames_to_column("Sample") %>%
as_tibble() %>%
gather(gene, cpm, -Sample) %>%
inner_join(tibble(
Sample = colnames(getData(cObjSng, "counts")),
Classification = getData(cObjSng, "classification")
)) %>%
mutate(type = case_when(
str_detect(Classification, "^S") ~ "Small intestine",
str_detect(Classification, "^C") ~ "Colon",
TRUE ~ "Other"
)) %>%
group_by(gene, Classification) %>%
summarize(mean = mean(cpm), pct = 100 * (length(which(cpm != 0)) / n())) %>%
mutate(scaled.mean.exp = scale(mean)) %>%
ungroup() %>%
mutate(type = case_when(
str_detect(gene, "^Wnt") ~ "Wnt ligand",
gene %in% c("Lgr5", "Axin2", "Ascl2") ~ "Wnt target",
TRUE ~ "error"
)) %>%
mutate(Classification = parse_factor(Classification, levels = classOrder.MGA('total'))) %>%
mutate(gene = parse_factor(gene, levels = rev(gene.order))) %>%
mutate(type = parse_factor(type, levels = c("Wnt ligand", "Receptor", "Wnt target"))) %>%
ggplot() +
geom_point(
aes(Classification, gene, size = pct, colour = scaled.mean.exp)
) +
facet_grid(type ~ ., space = "free", scales = "free") +
scale_color_gradient(low = "white", high = "darkred") +
scale.func(range = c(0, 18), limits = c(scale.min, scale.max)) +
guides(
size = guide_legend(
title = "% expressed", title.position = "top", title.hjust = 0.5
),
colour = guide_colorbar(
title = "Scaled mean expression", title.position = "top",
title.hjust = 0.5, barwidth = 10
)
) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 90, hjust = 1),
axis.title = element_blank(),
legend.position = "top",
legend.justification = "center"
)
## Joining, by = "Sample"
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.stemness.dotplot.pdf',
device = cairo_pdf,
height = 160,
width = 250,
units = "mm"
)
Relative frequency of cell types in singlets vs. deconvoluted multiplets
singlets <- c(table(getData(cObjSng, "classification")))
singlets <- singlets / sum(singlets)
deconv <- colSums(adjustFractions(cObjSng, cObjMul, sObj))
deconv <- deconv[match(names(singlets), names(deconv))]
deconv <- deconv / sum(deconv)
if(!identical(names(singlets), names(deconv))) stop("name mismatch")
p <- tibble(
class = names(singlets),
singlet.freq = singlets,
multiplet.freq = deconv
) %>%
ggplot() +
geom_point(aes(singlet.freq, multiplet.freq, colour = class), size = 3) +
scale_colour_manual(values = palette('total')[order(names(palette('total')))]) +
xlim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
ylim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
geom_abline(slope = 1, intercept = 0, lty = 3, colour = "grey") +
labs(x = "Singlet relative frequency", y = "Multiplet relative frequency") +
guides(colour = guide_legend(title = "Cell Type")) +
theme_few()
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.sngMulRelFreq.pdf',
device = cairo_pdf,
height = 160,
width = 250,
units = "mm"
)
plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 0, classOrder = classOrder.MGA('total'),
classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85
)
## Joining, by = "class"
Detected duplicates - quadruplicates.
Only ERCC estimated cell number max 4.
Weight cutoff = 10.
# adj <- adjustFractions(cObjSng, cObjMul, sObj, binary = TRUE)
# samples <- rownames(adj)
# rs <- rowSums(adj)
# keep <- rs == 2 | rs == 3 | rs == 4
plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 10,
classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
alpha = 1e-3
)
## Joining, by = "class"
pdf('../figures/MGA.enge.circos.pdf', width = 12, height = 10, onefile = FALSE)
plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 10,
classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
alpha = 1e-3
)
## Joining, by = "class"
invisible(dev.off())